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ABSTRACT 

We present the first three-dimensional, full-star simulations of convection in a 
white dwarf preceding a Type la supernova, specifically the last few hours before 
ignition. For these long-time calculations we use our low Mach number hydro- 
dynamics code, MAESTRO, which we have further developed to treat spherical 
stars centered in a three-dimensional Cartesian geometry. The main change re- 
quired is a procedure to map the one-dimensional radial base state to and from 
the Cartesian grid. Our models recover the dipole structure of the flow seen 
in previous calculations, but our long-time integration shows that the orienta- 
tion of the dipole changes with time. Furthermore, we show the development of 
gravity waves in the outer, stable portion of the star. Finally, we evolve several 
calculations to the point of ignition and discuss the range of ignition radii. 

Subject headings: supernovae: general — white dwarfs — hydrodynamics - 
nuclear reactions, nucleosynthesis, abundances — convection — methods: nu- 
merical 



1. Introduction 

Modeling highly subsonic convection in stars requires algorithms designed for long time 
integration. In the low Mach number approximation, we filter out sound waves while keep- 
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2006aJ — henceforth paper I, lAlmgren et al.ll2006bl — henceforth paper II, and 



2008| henceforth paper III), we developed a low Mach number stellar hy- 



drodynamics algorithm for reacting full-star flows in order to study the convective phase of 
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Type la supernovae (SNe la). In paper I, we derived the low Mach number equation set. 
In paper II, we included the effects of heat release due to external sources and allowed for 
a time-dependent background state. In paper III, we incorporated reactions into the system 
and also allowed the background state to evolve in response to large-scale convection and 
large-scale heating. Here, we extend the algorithm to spherical full-star problems using a 
three-dimensional Cartesian grid geometry. 

Our target application for this algorithm is the period of convection that precedes the 
ignition of SNe la. The standard model of a SN la involves a white dwarf in a binary 
system accre ting from a normal companion , and approaching the Chandrasekhar mass (see 
for example iHillebrandt fc Niemeyerl |2000| ). The increase in the central temperature and 
density accompanying the accretion seed carbon burning in the core, which in turn drives 
convection in the star. This convective 'si mmering' phase can last centuries, slowly increa sing 
the core temperature of the white dwarf (jWoosley et al.ll2004l ; IWunsch fc Woosleyl 120041 ) . A 
similar starting condition might be achieved in merg ing white dwarfs i f mass is added slowly 
enough to avoid ignition at the edge of the stars (jYoon et al.l 120071 ). During this phase, 
fluid heated by reactions buoyantly rises and cools via expansion, exchanging heat with 
its surroundings. The extent of the convective region grows with increasing temperature, 
eventually covering roughly the inner solar mass of the star. Outside of the convective region, 
the star is stably stratified. 

The continued increase in central temperature, coupled with the extreme temperature 
sensitivity of the carbon reactions, means that eventually the reactions proceed vigorously 
enough that a hot bubble cannot cool fast enough, and a burning front is born . This happens 
for a temperature of about 7 — 8 x 10 8 K ( jNomoto et al.lll984l ; IWoosleyill990l ). This burning 
front will quickly propagate through the white dwarf, converting most of the carbon/oxygen 
fuel to heavier elements, and releasing enough energy to unbind the star. However, ex- 
actly where in the star the ignition takes place is still unk nown. Among the earlie s t wor k 
to consider the role o f buoy ancy i n off-center ignit i on we re iGarcia-Senz fc Woosleyl (119951 ). 
Bychkov fc Libermanl (Il995l ). and iNiemeyer et al.l (119961 ). Additionally, some multidimen- 
sional studies have be e n done of the dynamic s of the first bubbles to ignite in a white dwarf 
( llapichino et al.ll2006l ; IZingale fc Dursil 120071 ). These papers made the case that we really 
need to understand whether the ignition is at the center or off-center. As calculations have 
become more sophisticated, it has only become more clear that the outc ome of the explosion 



is extremely sensit i ve to exactly how the burning fronts are initi ated (IGamezo et al.l 12005 
Jordan et al1l2008l : iRopke et al1l2007l : IGarcia-Senz & Bravd 120051 ) . 



Less work has been done on multidimensional modeling of the convective phase preceding 
the explosion. To date, no multidimensional calculation of the convection in the white 
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dwarf has modeled the entire star. The major contributions thus far are two-dimensional 



simulations of a 90° wedge of the star using an implicit hydrodynamics code (Hoflich fc Stein 



2002 



Stein fc Wheeler! 120061 ). and a three-dimensional anelastic calculation (IKuhlen et al 



20061 ) of the inner convective region of the star. All of these calculations cut out a small 
part of the central region of the star to avoid the coordinate singularity at the origin in 
spherical coordinates. Furthermore, the Kuhlen et al. calculation modeled the star out to a 
radius of only 500 km — leaving out p art of the convec t ive zo ne and the surrounding, stably 
stratified region. The calculations by IHoflich fc Stein! (120021 ) found ignition near the center 
of the white dwarf, produced by the fluid flow converging toward the center, with convective 
velocities of about 100 km s _1 . However, the ignition they see was likely affected by the 
converging ge o metry of their computational domain. The three-dimensional calculations by 
Kuhlen et al.l (120061 ) showed that the large-scale flow took on a dipole pattern, suggesting 
that off-center ignition in an outflow on one side of the star might be favored. Th e y also 
investigated the role of rotation. Finally, recent calculations shown in IWoosley et al.l (120071 ) 
used an anelastic method on a Cartesian grid, avoiding the singularity in the center, but still 
cut out the outer part of the convective region and the convectively stable region surrounding 
it. Here the dipole was once again seen. 

As seen from the wide range of explosion outcomes in the literature, realistic initial 
conditions are a critical part of SNe la modeling. Only simulations of this convective phase 
can yield the number, size, and distribution of the initial hot spots that seed the flame. 
Additionally, the initial turbulent velocities in the star are at least as large as the laminar 
flame speed ( IHoflich fc Stein! |2002| ). so accurately representing this initial flow may be an 
important component to explosion models. Perhaps o wing to a limited number of convection 
calculations, with few exceptions (ILivne et al.ll2005l ). nearly all explosion models to date 
begin with a quiet (zero velocity) white dwarf. 

Our goal in this study is to demonstrate that we have developed low Mach number 
hydrodynamics to the point where we can perform detailed calculations of the convective 
flow preceding the explosion, and to begin to understand the nature of the dynamics. In 
this work, we model the entire st ar, including the reg ion surrounding the convective zone. 
Recently, it has been suggested (IPiro fc Chang 120081 ) that the dynamics at the interface 
between the convective and stably-stratified regions of the star may be important during 
the flame propagation phase. Only full star calculations can capture this part of the flow. 
The resulting simulations can then form the basis for simulations of the flame propagation 
to build a more detailed picture of SNe la. 
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2. Numerical Methodology and Setup 

The basic idea of low Mach number hydrodynamics is to reformulate the fluid equa- 
tions to filter out sound waves while retaining the compressibility effects important to the 
problem — in this case, local compressibility effects due to burning, and large scale effects 
due to the background stratification of the star. A full derivation of the equations of low 
Mach number hydrodynamics is presented in papers I — III. Here we show the final equations 
and discuss adjustments needed for the spherical star. We recall that the use of low Mach 
number equations rather than the fully compressible equations enables the use of a time step 
based on the fluid velocity rather than the sound speed; this allows a 1/M increase in the 
time step over traditional compressible codes, where the Mach number, M, represents the 
ratio of fluid velocity to sound speed. During the convective phase preceding the first flames 
in SNe la, we expect the Mach number to be 0(0.01), making a low Mach number algorithm 
an appropriate choice. 

We choose to discretize our three-dimensional grid using Cartesian rather than spherical 
coordinates in order to avoid a coordinate singularity at the center of the star. This gives 
rise to the most notable difference from paper III — the base state is a 1-d radial profile and 
is not aligned with the any of the axes in the three-dimensional Cartesian grid. We refer 
to this as a spherical geometry, reflecting the fact that the base state is discretized in 1-d 
spherical coordinates. Throughout this paper, we refer to the Cartesian coordinates of the 
center of the star as (x c , y c , z c ). 



2.1. Equation Set 

The formulation of our equations relies on the existence of a base state density, Po( r )> 
and pressure, Po( r ) 5 that are in hydrostatic equilibrium, Vpo = Poge r , where e r is a unit 
vector pointing in the radial direction from the center of the star. In spherical geometries, 
the gravitation acceleration, g(r), is computed solely using the base state density as 

GM cncl (r) 

g(r) = - 2 — (1) 

with the mass enclosed within a radius r defined as 

M encl (r) = 4vr f p (r')r' 2 dr' . (2) 



As we discuss in § 12.41 we use a cutoff density, p cuto s, m our initial model. The star is mapped 
onto the grid down to this cutoff density, surrounded by an ambient medium. In computing 
Mend, we stop contributing to M enc i once the density drops below p cu toff- 
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In this paper, we reuse much of the notation from paper III. An overbar represents the 
average of a quantity over a layer of constant radius in the star, 



1 



A(n H ) 



(x) dA , (3) 



H 



where Qh is a region at constant radius in the star, and A(Qh) = Jq h dA. In this notation, x 
represents the Cartesian coordinates on the 3-d grid and r is the base state radial coordinate 
centered at (x c , y c , z c ). A subscript '0' represents a base state quantity. We compute e r in a 
cell indexed by k) with Cartesian coordinates (xi,yj,Zk) as 

%i X c Uj Vc z k / 4 % 

G r G x T" Gy T" G z , y±j 

ry ry ** ry 

with r 2 = (xi — x c ) 2 + (yj — y c ) 2 + (z k — z c ) 2 , and g x , g v , g z the unit vectors for the Cartesian 
coordinate system. 

In our previous work, the total fluid velocity, U, was decomposed into a local velocity 
field, U and base state velocity, w , as 

U = U(x,t) + w (r,t)G r . (5) 

The base state velocity, is used to adjust the base state in response to the heating on the 
grid. In paper II, we demonstrated that when the heating is large, expanding the base state 
is critical to accurately modeling the flow. 

In the current application, convection in the white dwarf, the heating is small until the 
flame ignites. Therefore, for these first calculations, we use a background state that is fixed 
in time. We will later quantify the extent to which this assumption of a fixed background 
state is valid. This simplifies the evolution equations, and we can now use U for U and 
w = 0. 



The full state evolves according to 
d( P X k ) 



at 



V-(UpX fe ) + P Cj k , (6) 



^ = -U.VU-iw-^},e r . (7) 
at p p 

Equation is the species evolution equation, where X k is the mass fraction of species k, with 
creation rate u k provided by the nuclear reaction network. The mass density, p, is simply 
p = ^2k(pXk)- For the velocity evolution equation, ©, n is the dynamic pressure resulting 
from the asymptotic expansion of the pressure in terms of Mach number. In paper III we 
also evolved the enthalpy, for the sole purpose of getting the temperature to feed into the 
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reaction network. For this paper, we instead define the temperature from p, po, and X k . 
Our experience has shown that, with the spherical geometry, the discretization errors are 
minimized by using the hydrostatic, radial base state pressure to define temperature. We 
will revisit this in a future paper. This system of equations is identical to that presented in 
paper III, with U = U, wq = 0, and dpo/dt = 0. 

The velocity field is subject to a constraint equation, 

V-(A)U) = foS , (8) 

with 

AW=M 0)exp(^f^<) , (9) 
where I\ is the average over a layer of d(logp)/d(log p) at constant entropy, and 

S = -a £ k U k + — Px k UJk + 0"#mic • (10) 

Here, p Xk = dp/dX k \ p>TXj ^ k , £ k = dh/dX k \ pTX ..^ k , p p = dp/dp\ TXk , and a = p T /(pc P p p ), 
with p T = dp/dT\ pXk and c p = dhj 'dT\ Xk the specific heat at constant pressure. In these 
derivatives, h is the specific enthalpy, defined in terms of the specific internal energy, e, 
pressure, and density as h = e +p/p. Finally, H nnc is the nuclear energy release (with units 
of erg g _1 s _1 ) as computed from our reaction network. Physically, S represents the local 
compressibility effects due to heat release from reactions and composition changes. The 
presence of the density-like quantity Po inside the divergence in the constraint captures the 
expansion of a parcel of fluid as it rises in the hydrostatically stratified star. 

We refer the reader to the extensive comparisons with compressible algorithms in papers 
I through III that demonstrate the validity of the low Mach number approximation. For the 
most part, the algorithm to evolve the star follows closely that described in paper III. For 
the construction of the advective terms, the interface states are again co nstructed using a 



piecewise linear unsplit Godunov scheme ba sed on tha t of IColellal (Il990l ). but we now use 



the full corner-coupling scheme developed by ISaltzmanl fjl994l ). In the subsections below we 



point out the differences for the present application. 



2.2. Mapping 

Since the one-dimensional radial base state is not aligned with any of the axes in the 
three-dimensional Cartesian grid, the discretization of quantities that involve both the base 
state and the full state becomes complicated. Various parts of the algorithm (such as the 
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averaging operations) require a mapping between the base state and the full state. Because 
the base state is not aligned with the Cartesian coordinate axes, we are free to choose the 
base state resolution independent of the Cartesian grid spacing. Numerical experimentation 
has shown that setting the base state resolution, Ar, to be finer than the Cartesian grid 
resolution, Ax, gives the best results. (Here we assume Ax = Ay = Az.) For the present 
simulations, we use 5Ar = Ax. We refer to the procedure that maps data from 1-d to 3-d 
as f ill_3d, and the complementary procedure that maps from 3-d to 1-d as average. 

Figure Q] shows the Cartesian grid overlaid by the spherical base state (for simplicity, 
the figure is drawn in 2-d using 2Ar = Ax). The f ill_3d procedure computes the distance 
of the center of cell indexed by (i, j, k) from the center of the star, 



Xi ~ x c ) 2 + (yj ~ Vc) 2 + (z k - z c ) 2 . (11) 

We use this radius to find the corresponding radial bin as n = int(r/Ar) (here, our conven- 
tion is to use 0-based indexing for the base state). We can then initialize a Cartesian cell 
quantity q from its corresponding base state quantity, q , as g^fc = qo, n - 

For the average process, we first define a coarse 1-d radial array with Ar c = Ax. 
Then, for each cell indexed by (i, j, k) we again compute the radius, r, as above, and define 
the index of the corresponding coarse radial bin, n c = int(r/Ar c ). We define go,n c as the 
average of all the qij,k whose Cartesian cell centers map into coarse radial bin n c . Next, we 
construct edge-centered states on the coarse radial bin using the fourth order approximation, 
Qo,n c +y 2 = (7/12)(<? ,n c + <?o,n c +i) - (V 12 ) (?o,n c -i + <?o,n c +2)- Finally, for each coarse radial 
bin we construct a quadratic profile using go,n c -y 2) <?o,n c and <Zo,ra c +y 2 - This is based on the 



interp olating; polynomial used by the PPM scheme to find edge states (IColella fc Woodward 



19841 ) . Specifically, for n c Ar c < r < (n c + l)Ar c , the interpolating polynomial is 

Qo(r) = <?o,n c -y 2 + £(r) {Aq nc + g 6 , n Jl - £(r)]} , (12) 

with 

r — n c Ar c , . 

t{r) = c c , 13 
Ar c 

A<?n c = <?0,n c +y 2 - <?0,n c -y 2 , (14) 

and 



<?6,n c = 6 



<?0,n c - ^ (<70,n c +y 2 + %,n c -y 2 ) 



(15) 



We note that since we are not evolving the base state in the simulation presented here, the 
feedback from the full state to the base state through average is limited to computing Ti, 
as needed for updating j3 . 
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2.3. Microphysics 



We use t he ge neral stellar equation of state described by iTimmes fc Swestyl (120001 ); 
Fryxell et al.l (120001 ). which includes contributions from electrons, ions, and radiation. For 
these calculations we include the effects of Coulomb corrections included in the publicly 
available version of this EOS (ITimmes! 120081 ). 



Our reaction network is u nchanged from paper I I I, and is a single-step 12 C+ 12 C reaction 
using s creening as describ ed in lGraboske et al.l (119731 ) ; IWeaver et al.l (119781 ) ; lAlastuey &; Jancovici 



( 119781 ): lltoh et all (119791 ). resulting in 24 Mg ash. We release the energy corresponding to the 
binding energy difference between the magnesium ash and carbon fuel. Paper III provides 
full details on how the reaction network is solved. Our only change from the implementation 
there is that we now update the temperature at the end of the reaction step. Finally we 
note that we do not call the reaction network for densities below p cu toff- 

We note t hat by integrat i ng th e reaction rate equation, we are dealing with reactions 
differently than lKuhlen et al.l (120061 ). There, an analytic approximation to the reaction rate 
was used and evaluated given a temperature and density. Our method extend s more easily 
to a full reaction network. A second difference is that iKuhlen et al.l (120061 ) burned to a 
mix of neon and magnesium, leading to a slightly lower energy release. This difference may 
affect the timescales we see in the calculation, but we don't expect it to introduce qualitative 
differences. 



2.4. Initial Model 



We begin with an initial 1-d white dwarf model produced with the stellar evolution 
code, Kepler (IWeaver et al.lll978l ). This model was evolved to the point where the central 
temperature is 6 x 10 9 K and the central density is 2.6 x 10 9 g cm -3 . The composition is 
about half 12 C and half 16 0, with a small amount (< 0.5%) of ash in the center of the star. 
The total mass of the star is 1.382 M Q . 



We follow the procedure outlined in IZingale et al.l (120021 ) to convert the initial model 
from the one-dimensional Lagrangian mesh used by Kepler to the uniformly-zoned Eulerian 
grid used in our calculation. It is important that the initial model satisfy hydrostatic equi- 
librium discretely with our equation of state on the base state grid we use for our simulation. 
In particular, we want to enforce the following discretization of hydrostatic equilibrium: 

1 



PO,: 



PO,: 



-Ar(p ,i + po,i+i)9i 



+y 2 



(16) 



Hydrostatic equilibrium alone does not specify our initial model, since we must also specify 
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the initial temperature. In the interior of the star, where convection dominates, constant 
entropy is a good approximation. We use this constraint together with equation ffl6|) and 
the equation of state to find the temperature, density, and pressure throughout the inner 
region of the star. For the composition, we use the profile provided by the Kepler model, 
but, since we are using a reduced network, we group together the 20 Ne and 24 Mg ash into a 
single composition variable. 

The convective region is surrounded by an outer, convectively stable region. When the 
isentropic temperature profile drops below the temperature provided by the Kepler model, 
we switch to using the Kepler temperature. Figure [2] shows our final temperature profile, 
along with a completely isentropic model for reference. The departure of the two temperature 
curves marks the boundary of the convective region. The mass of the inner isentropic region 
of the star is 1.131 M . We note that the spatial extent of the convective zone in the white 
dwarf is somewhat uncertain. Different assumptions about the accretion history of the white 
dwarf would lead to different mass convection zones. 

Overall, this procedure results in a slight adjustment of the structure of the star com- 
pared to the initial Kepler model. The resulting model serves as the initial base state for our 
calculation. As discussed in papers II and III, outside of the star we cannot bring the density 
down to arbitrarily small values, as that would result in too high a velocity (a consequence 
of our constraint equation). In practice, we impose a cutoff at a moderately small density, 
Pcutoflh and set the density to this constant value outside of the star. For the main calculation 
presented here, we choose p cu toff = 3 x 10 6 g cm -3 . While this may sound high, we note that 
the mass of the star enclosed by p cu tofr is 1.378 M & — a 0.2% difference from the total mass 
of the star. We note also that as in paper III, we use an anelastic cutoff, the density below 
which the coefficient, /? , of our velocity constraint is defined by keeping Po/po constant. In 
this paper, we always set the anelastic cutoff to be p C utoff- 

The initial three-dimensional state is set by using the fill_3d routine in § 12.21 to inter- 
polate po, po, Xkfi, and T to each cell center. The initial velocity field is not as well defined. 
The one-dimensional stellar evolution code used mixing length theory to describe convective 
mixing in the interior of the star. When we map the model onto our three-dimensional grid, 
there is a region that is convectively unstable (corresponding to the region in Figure [2] where 
r < 1.0 x 10 8 cm). However, there is not enough information in the one-dimensional model 
to initialize a three-dimensional velocity field that correctly represents the convective field. 

If we start with zero initial velocity, then at t = the reactions near the core generate 
a large amount of energy and the highly nonlinear form of the reaction rate means that the 
energy release quickly grows. Without an initial velocity field to advect some of this energy 
away from the core, the energy generation grows too quickly, and an unphysical runaway 
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occurs. However, by starting with an initial non-zero velocity field, our simulation very 
quickly finds a convective velocity field that balances the energy generation at the core. 
Thus we define a set of Fourier modes, 



l,m,n 



r (y) 

l,m,n 
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l,m,n 
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sin 
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a 

( 2imz 



Lm,n 



(17a) 
(17b) 
(17c) 

(18a) 
(18b) 
(18c) 



where a is the characteristic scale of the perturbation and the 4>i r ^^ are randomly generated 
phases between [0, 2tt]. We then compute the total contribution to the velocity perturbation 
from the modes as 

3 3 3 



u 



1=1 m=l n=l N hm,n 
3 3 3 



W 



1=1 m=l n=l 

3 3 3 

ZEE 



1 



1 ^> m >™ 



Hl,rn,n h ^ i m ,n Lm,n Lm,n ' u l,m,n' n ^Lm,n^Lm,n J Lm,n 



where ai irnin , A,m,n> and 7i, m ,n are randomly generated amplitudes between [-1, 1], and 
-^,m,n = a// 2 + m 2 + n 2 is the normalization, 

A perturbational velocity field is then computed as 

Au' ~ 



u 



2 
Av' 
~Y 
Aw' 



l + tanh( rpert r 
a 

l + tanh( rpcrt ~ r 



1 + tanh 



'"pert f 
d 



(20a) 
(20b) 
(20c) 
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where the tanh profile gradually cuts off the perturbation at a radius r pert with a transition 
thickness d. Finally, the initial velocity field is computed by applying the projection to 
(u", v", w") to ensure that it satisfies the divergence constraint. We pick the amplitude, A, 
to be small, and independent of the velocities used in the 1-d stellar evolution model. Once 
the flow field is established, we expect the details of the initial velocity field to be forgotten. 
This is an area we will explore in a subsequent paper. 

Throughout the calculation we solve the reaction network to compute the energy release 
that drives the convection. By starting at a low initial central temperature, we thus expect 
a realistic flow field to build up over time as the central temperature in creases from the reac - 
tions. In this respect we differ from the initialization procedure used in lKuhlen et al.l ( 120061 ). 
In their anelastic approximation, they carried the perturbational temperature separately 
from the base state temperature, and to initialize the flow field they evaluated the carbon 
burning heating term using only the base state temperature. By excluding the perturbational 
temperature, they left out the nonlinear feedback in the extremely temperature-sensitive car- 
bon reaction rate, and therefore built a flow field without the chance of runaway. Once the 
flow field was established, they fed the temperature perturbations back into the reaction rate 
to watch the runaway. 



2.5. Sponging 

As described in paper III, we use a sponge to damp the velocities outside of our region 
of interest. We use the same functional form here, with the velocity forcing given by 



U old - A^/ damp U r 



where k is a frequency. For all results presented here, we use k 
has the form: 



(21) 

10 s -1 . The sponge factor 



damp 



cos 



7T 



r — r. 



sp 



tp 



sp 



if r < r 



sp 



if r sp < r < r tp 
if r > r tp 



(22) 



The quantity r sp represents the radius where the sponging term gradually begins to turn on, 
and is set to the radius corresponding to 10-p cuto fr- The top of the sponge, r tp , where the 
sponge is in full effect is set as r tp = 2r m d — r sp , with r m( j set to the radius corresponding 
to the Pcutoff- As noted above, we use a p C utoff = 3 x 10 6 g cm -3 for these calculations, 
so the corresponding density where our sponging begins is 3 x 10 7 g cm -3 . Based on our 
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initial model, 1.320 M of the star is contained within r sp — the sponge only affects the 
very outer portion of the star. Figure [2] shows the location of r sp for our initial model — we 
see that it is well outside the convectively unstable region. Figure [3] shows /damp vs. r for 
Pcutoff = 3 x 10 6 g cm" 3 . 

This sponge is effective in damping the velocities at the edge of the star. Our domain 
is D = 5 x 10 8 cm on a side, so the distance from the center of the star along one of the 
coordinate axes to the edge of the domain is 2.5 x 10 8 cm. The distance from the center to a 
corner of the domain is v3 larger. Because we are placing a spherical star in a cubic domain, 
we found that we need an additional sponge to damp the velocities in the outer corners of 
the domain — well outside of the star. We define an outer sponge of the same form as above, 
but with r tp = D/2 and r sp = r tp — 4 Ax, where Ax is the grid spacing, and k set to 10 times 
the value of the inner sponge. This additional sponge is included in the momentum equation 
in the same fashion as the inner sponge. Figure [3] shows the profile of this additional sponge 
as well. 



3. Results 

Our main goal in these simulations is to study the convection in the white dwarf up 
to the point of ignition. In this section we present results for our main 384 3 convection 
calculation, supporting calculations with lower resolution, as well as a test problem. In each 
case, the code was run with an advective CFL number of 0.5 with the star centered in a 
domain 5 x 10 8 cm on a side. 



3.1. Test Problem: Isentropically Stratified Star 

To test the interaction between the spherical base state and the 3-d Cartesian repre- 
sentation of the star, we perform a simple advection test with an analytic solution. First, 
we construct a completely isentropic initial model. This is achieved by picking a central 
density of 2.6 x 10 9 g cm -3 and a central temperature of 6 x 10 8 K, and a uniform com- 
position of 0.3 12 C and 0.7 16 0, and integrating outward using our hydrostatic equilibrium 
constraint, equation ffT6l) . and forcing the entropy to be constant through the equation of 
state. We initialize the full state using the isentropic base state with no perturbations. We 
also set /3 = p discretely (which is true analytically for an isentropic base state and constant 
Ti), and disable all reactions and heating. The constraint is now identical to the anelastic 
constraint, V • (poU) = 0. 
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Under these conditions, the continuity equation becomes: 



dp 
di 



V ■ (pU) 



V • ( Po U) = 



(23) 



using the fact that p = po initially, and the anelastic constraint. As a result, we see that the 
density should remain constant in the star regardless of the velocity field. 

This provides a means to test our mapping procedure. If we start with an isentropically 
stratified star and seed a random velocity field, the density should not change with time. 
For our test, we start with a random velocity field described by equation fl20l) . For the 
amplitude of the perturbation, we set A = 10 7 cm s _1 — this is typical of the highest velocities 
we expect to see in our convection calculations. For the size of the perturbation, we set 
r pert = 5 x 10 7 cm — this value represents about half the size of the expected convective 
region in the white dwarf. Finally, we set the characteristic wavelength of the perturbation, 
a = 10 7 cm. We make the transition between the perturbation and the ambient star sharp, 
effectively smaller than our grid resolution, setting d = 10 cm. The resolution is 384 3 , the 
same as that used in the main calculation in the next section. 

To assess the change in density with time, we will look at the average density, (p), as 
a function of radius, and the deviation of the density as a function of radius, dp. We define 
these as 



where Q r is the set of cells in the computational domain whose center falls within the radial 
bin at radius r, and Nq t is the number of cells in Q r . The RMS fluctuations are 



Here we recognize that the base state density, po represents the average density at a given 
radius. We compute and store (p — po) for every zone in our computational domain directly 
in the code as the simulation runs, and then compute (Sp) r using equation ( |25i) with a radial 
bin spacing Ar = Ax — this ensures that no interpolation is needed to fill radial cells. 

Figure H] shows a plot of {Sp) r / (p) r vs. r at several times. By normalizing to the average 
density, (p) r , we are seeing a measure of the relative error in the density from our advection 
scheme. As the plot shows, even after 1500 s of evolution, the error at the center of the 
star is < 10~ n . Once we are outside of the star r > 2 x 10 8 cm, the error rises, but still 
stays below 5 x 10~ 9 everywhere. This demonstrates that our algorithm accurately preserves 
dp/dt = in the limiting case of an isentropic model with no heating and /?o = Po- 




(24) 




(25) 
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3.2. Convection in a White Dwarf 

We model convection in the white dwarf by mapping the initial model described in § 12.41 
onto our Cartesian grid. For the initial velocity field, we use A = 10 5 cm s _1 , d = 10 5 cm, 
r P crt = 2 x 10 7 cm, and a = 10 7 cm. 



3.2.1. Diagnostics 

To help us understand the character of the flow in our calculations we make use of 
several diagnostic quantities. We define the region of interest of the domain, fistar, to be 
those computational cells with p > p cu tofr, where we have used p cu tofr = 3 x 10 6 g cm -3 unless 
otherwise specified. The diagnostics defined below are computed every time step, as the 
code is running. 

The peak temperature in the domain is simply 

^peak = max{T} . (26) 

i 'star 

As the temperature in the star increases considerably toward the center of the star, we expect 
the peak temperature to be close to (but not exactly equal to) the central temperature. 



Motivated by previous results that suggest a dipole nature to the flow (IKuhlen et al. 



20061 ). we look at several diagnostics based on the radial fluid velocity. First we define the 
radial velocity to be v r = U-e r . Then we compute components of the density-weighted 
average radial velocity in each coordinate direction, 

(v r )x = Yl pVr { r J / P ' ( 27 ) 

where r is the distance of a given zone from the center of the star, and iVn atar is the number 
of computational zones contained in the domain O s tar- We similarly compute (v r ) y and (v r ) z 
using the y and z coordinates and centers. The relative magnitudes of the components of 
(v r )i tell us about the direction of any dipole-nature to the flow. In particular, we can derive 
the directional angles in the x-y plane, and 9 as measured from the z-axis as 

\\ v r/x/ 

and 

9 = tan" 1 I -5! -— I . (29) 
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We could have instead computed the average radial velocity without a density weighting, 
but because we are summing over the entire star (where p > p C utoff) ; we are including the 
outer convectively stable region in the average, where we do not expect to see much influence 
from the dipole. By density-weighting, we are giving more weight to the center of the star, 
where the convective pattern dominates. 

Finally, to get a sense of scale for the radial velocity in the star, we compute 

Or) P eak = max{|tv|} . (30) 



3.2.2. Long Term Convective Behavior 

Our main result is a 384 3 calculation of convection in a white dwarf, starting from an 
initial model with a central temperature of 6 x 10 8 K. Our computational domain is 5 x 10 s cm 
on a side, giving us 13 km zones. Our goal is to follow the convection as reactions bring the 
central temperature up over 7 x 10 s K, and into the regime of ignition. 

As noted in § 12.41 we started with a small velocity perturbation near the center of the 
star and the velocity otherwise zero. As the simulation begins, reactions heat the core of the 
star, and since the background of the star is isentropic, the heated fluid at the core begins to 
buoyantly move radially outward. Figure \5\ shows the magnitude of the vorticity (|V x U|) 
in the three orthogonal slice planes through the center of the star at several different times. 
At early times, we see convective flow developing near the center of the star. By 400 s, 
the convective flow has grown to fill the convectively stable region, and we see gravity waves 
excited in the stable region above. The later times show the convective pattern continuing to 
strengthen, with small asymmetries in the the vorticity moving through the inner convective 
region. For most of the simulation, we see a sharp distinction in the character of the flow 
at the boundary of the isentropic region in the star. However, toward the very end of the 
calculation, as shown in the very last pane of Figure we no longer see the separation 
between the two regions, and the convective plumes appear to travel through the entire star. 

Figure [6] shows contours of the radial velocity at four different times. Qualitatively, these 
times represent the early period (panel a, 800 s), two intermediate snapshots (panels b and 
c, 3200 and 3420 s respectively), and the very late stage of the calculation (panel d, 7132 s). 
Red indicates fluid moving radially outward, and blue indicates fluid moving radially inward. 
The gray surface is drawn at a constant density (p = p cu toff) an d represents the surface of 
the star. Very early we see the distinct asymmetric nature to the fl ow characteris t ic of a 



dipole flow. The dipole is not nearly as symmetric as that shown in iKuhlen et al.l (120061 ) 



perhaps due to differing resolution or the inclusion of the stably-stratified layer surrounding 
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the convective region in our study. In general, the outward moving fluid appears more 
coherent then the inward moving fluid. Comparing the images at different times, we see that 
the dipole direction changes with time. Occasionally, the flow takes on a more organized 
form, with the inward moving fluid forming a concentric ring around the outward flow, as 
shown in panel c. At the very late stages of the simulation (panel d) , we see what appears 
to be a breakdown in the distinction between the stable and unstable regions, with the flow 
much less organized and filling most of the volume of the star. The narrow gap between the 
velocity contours and the surface of the star at late times arises from our sponging term. We 
look at the sensitivity of the results to the position of the sponge in the next subsection. 

To get a better feel for the change in direction of the dipole, we compute the spherical 
angles, 9 and 4>, from (v r )i, as defined above. Figure [7] shows these angles as a function of 
time. We see that both angles move through their full range many times over the course of 
the simulation. We see that the characteristic timescale for to complete a circuit through 
27r is between 500 and 1000 s. At late times, it appears that the dipole is changing direction 
with a faster period, especially in the 6 plot. 

Figure M shows the peak radial velocity, (iv) pea k, inside the star, as a function of time. 
We see that it slowly rises with time, with a typical peak radial velocity of ~ 10 7 cm s _1 . 
Taking the convective region to have a radius of -R C onv ~ 10 8 cm, we define a lower bound to 
the convective turnover time of 2i? conv /(tv) pca k = 20 s. 

Ignition will occur when the reactions proceed so strongly that hot, reacting bubbles are 
not quenched by adiabatic expansion in the convective motions carrying the fluid away from 
the center of the star. Since the 12 C + 12 C reaction rate is so strongly temperature-sensitive, 
the peak temperature in the star serves as a good guide for observing the progression toward 
ignition. Figure M shows the peak temperature as a function of time for this calculation. We 
see a short transient at the start of the calculation where the temperature quickly rises and 
then settles back down — this occurs from the nonlinear feedback of the temperature into the 
reactions when the flow field is not yet fully developed. After a short amount of time, a 
convective flow field develops that properly matches the energy generation at the center of 
the star, and the temperature settles into a long, gradual rise. About halfway through the 
calculation, we can clearly see that the temperature rise is non-linear, and the temperature 
increase accelerates toward the very end, up to the point of ignition. The inset in Figure [9J 
shows the behavior of T pea k during the last 200 s. 

The reactions dump energy into the star, and it heats up throughout. Figure HOI shows 
the average temperature at a given radius as a function of radius at several different times. 
As we see, the temperature increases throughout the convective region. At late times, we see 
a distinct change in the temperature structure at the boundary of the convective region. This 
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change in temperature structure corresponds to the penetration of the vortical flow through 
the original boundary between the stable and unstably-stratified regions in the vortical plot 
(Figure [5]) shown above. It is not clear how robust this change in the character of the flow is 
to resolution — that is something that will be explored through higher resolution studies in 
the future. It is also the case that those outer layers, near the transition to a stably-stratified 
fluid, are where we would expect the expansion of the star to be greatest, so we need to check 
if neglecting the base state evolution was warranted. Figure [TT1 shows (Sp) r /(p) r vs. radius 
at several times. This is a measure of how much expansion of the star has taken place. If 
(Sp) r is large compared to (p) r , then the full state is carrying the expansion rather than 
the base state, and as we've shown in paper II, this can lead to inaccuracies. As Figure [11] 
shows, (Sp) r / (p) r is always below 1%, indicating the departure from the base state is small, 
and any expansion would be minimal. In each case, the curve at 7132 s corresponds to the 
point is ignition, discussed below. 

We can also look at the total kinetic energy in the star, which we compute as 



At the point when the peak temperature reaches 8 x 10 8 K, the total kinetic energy inside 
the star is 6.24 x 10 46 erg. To put this in context, we can compare the gravitational potential 
energy of the star, defined from our base state as 



with dM = 47ir 2 p dr. For our model, the gravitational potential energy is —3.2 x 10 51 erg. 
The internal energy of the gas is also quite large, E int = 2.7 x 10 51 erg, giving an energy 
difference of ~ —5 x 10 50 erg that needs to be overcome to unbind the star. Therefore, kinetic 
energy release up to the point of ignition is a tiny fraction of what is needed to unbind the 
star — as expected. 

The 384 3 calculation took 113156 time steps to reach a simulation time of 7131.8 s — at 
which point the peak temperature had risen to 8 x 10 8 K, and ignition shortly followed. 
Overall, the average time step is 0.063 s. At this same instant, the Mach number, attained 
in the outer layers of the star, reached a value of 0.079. Earlier in the calculation, the 
maximum Mach number in the domain was considerably lower. For comparison, the highest 
sound speed in the star (at its center) is 9.5 x 10 8 cm s _1 , which would give a corresponding 
time step of 7 x 10~ 4 s (assuming a CFL number of 0.5, and \U\ <C c s , where c s is the sound 
speed). 




(31) 




(32) 
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3.2.3. Effect of p cuto s 

As we noted above, the behavior of the coefficient in our constraint term and the location 
of the sponge are set by the density we refer to as p C utoff- To assess the influence of our 
choice of p cu toff; we perform a pair of simulations on a 256 3 grid that are identical except 
for the value of p cu toff- For our control case, we use p cuto fr = 3 x 10 6 g cm -3 , the value 
chosen for our main calculation. To explore the effects of lowering p cu toff> we also try a 
value of Pcutofr = 10 6 g cm -3 . We note that the location of the sponge in the momentum 
equation remains keyed to the choice of p cu toff , so with the lower value of p cu toff, the location 
of the start of the sponge moves outward from the center of the star. In terms of mass, 
Pcutoff = 10 6 g cm -3 means that the mass of the star enclosed is 1.381 M , compared to 
1.378 M & with p C utoff = 3 x 10 6 g cm" 3 . The location of the start of the sponge contains 
1.363 M , compared with 1.320 M with p cu toff = 3 x 10 6 g cm -3 . 

Figure [T2] shows T pca k as a function of time for the two calculations. As we see, the 
two curves track very well, indicating that the choice of p C utoff has little influence on the 
temperature behavior near the center of the star. The time at which final ignition occurs 
differs between these two cases by only 38.4 s out of over 6000 s of evolution. 

As noted above, the choice of p cut off is used to prevent the velocities from growing too 
large as the fluid experiences the steep density gradient at the edge of the star. We note 
that the time step the code takes with p C utoff = 3 x 10 6 g cm -3 is 23% larger then with 
Pcutoff = 10 6 g cm -3 . Thus it is computational favorable to use the slightly higher value of 
the cutoff density. 

3.2.4- Ignition 

As the inset in Figure M shows, up to the point of ignition, the peak temperature rises 
rapidly, only to fall again, as a spark fails to ignite. Figure [13] shows the temperature 
structure in the last 500 s for the 384 3 calculation and both 256 3 calculations, shifted so 
the time of ignition lines up. All three runs show the peak temperature fluctuating rapidly 
before ignition, indicating some hot spots failed to ignite. In each case, eventually, a hot spot 
burns faster than it cools and the temperature rapidly shoots up to over 10 10 K — igniting 
the first flame. At this point, our algorithm cannot deal with the rapid energy release, and 
the low Mach number approximation breaks down, so we stop the calculation. Physically, 
at this point the nuclear burning timescale is much shorter than the advection timescale. 
Whether a second hot spot ignites shortly following this one is not something we can address 
presently. 
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At the point of ignition, the radial velocity (Figured]) rises rapidly, reaching unphysically 
high values post-ignition — a symptom of the breakdown of the low Mach number method 
when the first flame ignites. However, prior to ignition, the algorithm remains valid, and 
we see that the radial velocities rise to ~ few xlO 7 cm s -1 . These velocities will affect the 
dynamics of the first flames to ignite. 

To determine the location of the ignition, we compute the radius of the hot spot that 
ignited by looking at the peak temperature. We define the ignition radius as the position of 
the hot spot at the time when the peak temperature passes beyond 8 x 10 8 K. In all cases, 
once this temperature is exceeded, the temperature quickly shoots up to O(10 10 ) K. 

For the high- resolution calculation, we find the location of ignition to be 21.6 km from 
the center of the star. The location of the peak temperature remained steady from (and 
possibly before) about 0.8 s before the peak temperature rose above 8 x 10 8 K, indicating 
that the hot spot was not moving very fast. The radial velocity in the zone with the peak 
temperature at the time we satisfy our ignition criteria is only 4.8 km/s. We note that 
physical center of the star is on a vertex on our Cartesian grid, so central ignition in this 
simulation would be at \/3Ax/2 = 11.3 km, the distance from the center of the closest grid 
cell to the center of the star. 

Figure [TH shows the perturbational temperature (full state temperature, T, minus the 
average temperature at the corresponding radius, (T)) in the central 64 3 zones (833 km on a 
side). We've picked the location of the slice planes to cut right through the hot spot. As the 
figure shows, the peak temperature is strongly localized to a single zone, with an extended 
hot region surrounding this location. 

For the two medium resolution (256 3 ) cases, we find the location of the hot spot when 
T P eak crosses 8 x 10 8 K to be 84.5 km (p cu toff = 3 x 10 6 g cm -3 ) and 32.4 km (p C utofr = 
10 6 g cm -3 ). For reference, central ignition at this resolution would correspond to 16.9 km. 
For the p cu tofr = 3 x 10 6 g cm -3 case, the location of the peak temperature changes rapidly 
(moving outward from the center) as the peak temperature crosses 8 x 10 8 K, ranging from 
32.4 km when T pcak = 7.51 x 10 s K at 6267.6 s to 89.0 km when T pcak = 8.13 x 10 s K at 
6269.9 s. The radial velocity in the zone satisfying our ignition criteria is 39 km/s. Clearly, 
the flow dynamics at the location of ignition are significant, in this case. We do not see 
movement of this magnitude for the p C utofi = 10 6 g cm -3 case, where the radial velocity in 
the zone satisfying our ignition criteria is only 2.9 km/s. For future calculations we will 
store the location of the hot spot along with the peak temperature at every time step to help 
better understand these dynamics. 

Taken together, we see a distribution of ignition radii in our results, ranging from 
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21.6 km to 89.0 km. Ignition is a highly nonlinear process, and we expect that were we 
to perform more runs, slightly tweaking the initial conditions, we would observe different 
values, all of which sample the distribution function of possible ignition locations in the 
problem. Owing to the stochastic nature of the problem, to really understand the ignition 
process requires performing a large number of slightly different calculations to map out the 
distribution function. 



3.2.5. Effect of Resolution 



It has been suggested that the behavio r of connective flow can dramatically change in 
character in flows of high Rayleigh number (lKadanoffN200ll ). In our simulation code, we do 
not explicitly add viscosity to the momentum equation (eq. [7j), so our Rayleigh number is 
determined by the numerical viscosity inherent in our advection scheme. Furthermore, the 
nature of the turbulence will depend on the Reynolds number of the flow, which again in our 
simulations is determined by numerical viscosity. Practically speaking, the way to increase 
the effective Reynolds and Rayleigh numbers of the simulation is to move to higher-order 
advection methods and to increase the resolution. 

While no amount of resolution will bring our effective Reynolds and Rayleigh numbers 
up to t he O(T0 14 ) and O(1 25 ) values, respectively, we expect in the true convecting white 



dwarf (jWoosley et al.ll2004l ). it is interesting to look at how the general results change with 
resolution. A second reason to explore resolution is that it is not known what size region 
will ignite. One might imagine that a region the size of only a few flame thicknesses needs to 
heat up to i gnite a flame. At the centr al densities in the white dwarf, the flame thickness is 
O(10~ 4 cm) (ITimmes &: Woosleylll992l ) — this is far below any resolution that can be obtained 
by a large-scale simulation code. However, the flame will initially burn in place, growing 
until it is lar ge enough (about 1 km ) that buoyancy becomes significant and it begins to rise 
and deform (IZingale fc Dursill2007l ). While this is still a smaller length scale than considered 
here, it is not out of reach with mesh refinement and larger computers. 

To begin to understand the effect of resolution, we consider three cases: 128 3 , 256 3 , and 
384 3 , corresponding to physical zone sizes of 39.1, 19.5, and 13.0 km respectively. We note 
that for the present study, computer resources prevent us from considering a 512 3 or higher 
case. 

Figure [15] shows T pea k vs. time for the three different resolutions. Immediately we see 
that the 128 3 case reaches ignition much faster than the two higher resolution cases. In fact, 
from our starting temperature of 6 x 10 8 K, the highest resolution run takes more than twice 
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as long in simulation time to reach ignition. Also apparent in the coarsest resolution run is 
that the temperature did not drop after the initial transient — this contributes to the faster 
overall evolution. Both of the higher resolution cases see a drop in the temperature after 
the initial transient, as the developing velocity field carries the heated fluid away from the 
center of the star. Close to ignition, the 256 3 and 384 3 runs show a similar slope in the T pca k 
vs. t curve shown in Figure [131 The peak radial velocity as a function of time also shows 
differences between the resolutions, as shown in Figure [161 There does appear to be some 
convergence with resolution. 



4. Conclusions and Discussion 



We have demonstrated that our simulation code, MAESTRO, is capable of following 
the convective flow in a white dwarf leading up to the ignition of a Type la supernova. We 
have explored the sensitivity of the results to resolution and to the choice of low density 
cutoff, p C utofr- Our test problem shows that discretizing the star on a Cartesian grid with a 
radial base state leads to an accurate representation of the flow. 

Over many convective turnover times, our simulations capture the rise of the peak 
temperature in the white dwarf up to ig nition, recover the dipole nature of the convective 
flow first shown in iKuhlen et al.l (120061 ). and track the change in direction of the dipole. 
We see, for the first time in multidimensional simulations, the distinct chang e in the nature 
of the flow at the outer boundary of the convective region, as discussed in iPiro fc Chang 
( 120081 ). The late time breakdown of this interface needs further investigation. 



All of our models reached ignition. For the two medium resolution runs, the ignition 
occurred at a radius of 32.4 and 84.5 km. For the high resolution run, it occurred at 
21.6 km. These are the locations of the first flames. We note that this is a highly-nonlinear 
problem, and small changes in the state of the star could affect the ignition process. To really 
understand the statistical distribution of initial ignition points requires running a suite of 
calculations, varying the initial model (central density, size of initial convective region), and 
the initial state of the star. With an ensemble of such calculations, we could get a much 
better understanding of the ignition process. We also need to understand how the ignition 
process differs with higher resolution. 

A detailed comparison to IHoflich fc Stein (hoo2h or IKuhlen et al.l ( 120061 ) is difficult, 
because of the differing geometries used. IHoflich fc Stein J2002h (hereafter HS) simulated 
a 90° wedge in 2-d, but cut out the innermost 13.7 km (in their "extended computational 
domain" run). As we discussed, in the present calculation the burning is strongly peaked 
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near the center of the star, so cutting out the center would miss a great deal of the energy 
generation. It would also prevent the fluid from flowing through the center of the star, which 
is the dominant pattern seen in the present calculation. Because HS used a spherical grid, 
and allowed the radial spacing to vary, there is no single grid resolution for their simulation, 
but they state that near the inner boundary, the grid resolution is "~ 2 km" — about 6x 
finer than the uniform resolution we use throughout the star. In both their calculation and 
our calculation, several hours of the runaway are followed leading up to the point of ignition 
(~ 3 hours for HS, ~ 2 hours for our calculation). Also, in both cases, the ignition takes 
place in a single zone. However, the details of the ignition differ — because HS have an inner 
boundary and a wedge-shaped domain, compression is generated that leads to the ultimate 
ignition near the center. In our calculation, we have flow through the center throughout the 
simulation. HS found the ignition to take place at a radius of 27 km in their model, which 
is within the range we report in our study. Furthermore, because their initial model had 
a strong gradient in the carbon mass fraction, with the outer portion of the star having a 
carbon mass fraction of 0.4 and the core having a value close to 0.25 (see HS, figure 1), HS 
report that the expanding convective region results in an increase in the carbon mass fraction 
at the center from 0.25 to 0.36. This is not the case in our model, where the carbon burning 
caused a 0.5% decrease in the central carbon mass fraction over the time we modeled. HS 
quote convective velocities between 40 and 120 km s _1 at ignition. In our case, the radial 
velocities were around 100 km s _1 for most of the evolution, rising to several times that just 
before ignition. 



The calculation by lKuhlen et al.l (120061 ) (hereafter KWG) modeled the star in 3-d, with 
an inner boundary at 50 km, again cutting out the center of the star. They also put the 
outer boundary at 500 km, which is approximately halfway through the convectively unstable 
region in the initial model (KWG and the present study use very similar initial models, 
generated from the Kepler stellar evolution code). Despite cutting out the center, KWG 
found that a large-scale dipole flow pattern dominates — similar to what we see (although we 
note that they also did a rotating model, and saw the dipole break down). Unlike HS or the 
present calculation, KWG do not model the evolution continuously leading up to ignition, 
but rather model two snapshots in time, corresponding to central white dwarf temperatures 
of 7 x 10 8 K and 7. 5 x 10 8 K, with dura tions of 70 and 41 s respectively. This amounts to a 



few turnover times (IKuhlen et al.l 120061 1 . It is difficult to compare resolution with KWG, as 
they use a spectral method for the spatial discretization. KWG quote typical velocities of 50 
- 100 km s _1 — consistent with both HS and the present calculation. Finally, in contrast to 
both HS and our calculation, KWG did not follow the evolution to the ignition of the first 
flame, but rather inferred from the size of the dipole flow pattern that ignition would likely 
be off-center. Overall the comparison to these previous calculations, and the variation seen 
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in our own set of calculations, suggest that more full-star, three-dimensional calculations, 
with varying initial parameters, are needed, to fully understand the ignition process. 

Future work will include both algorithmic improvements and more realistic physics. 
The focus of the next set of calculations will be to continue exploring the nature of the 
convection. On the algorit hmic front, we will switch to an unsplit impleme ntation of the 
piecewise parabolic method f Colella fc Woodwardll984 : Miller &: Colella 2002 ) for the advec- 
tion scheme and begin to incorporate adaptive mesh refinement. Together with an increase 
in grid resolution, these changes will allow us to push to higher effective Reynolds numbers. 
Physically, we will improve the nuclear energetics, add an enthalpy equation to better define 
the temperature throughout t he simulation, incor porate the exp ansion of th e base state, and 
include rotation. Numerical (IKuhlen et al.l 120061 ) and analytic (jPirol 120081 ) work has shown 
that the effects of rotation can be significant. 

The present calculations show only the ignition of the first flame, but the convection con- 
tinues and it is likely that other flames will ignite. This process of ong oing ignition could be 
critic al to the understanding of SNe la. To date, only limited studies (ISchmidt &: Niemeyer 
20061 ) have been performed investigating the effects of temporally-spaced ignition spots. By 
capturing the flames that ignite in our simulations and propagating them in a controlled 
fashion, we can continue a convective calculation to simulate the formation of additional 
ignition spots. 

Presently, our algorithm is unable to follow the evolution past the point where igni- 
tion occurs. However, right up until that point, the model remains valid. Therefore, these 
models can still provide useful starting conditions for explosion models run with fully com- 
pressible codes, simply by mapping the fully convective state right before ignition into a 
3-d compressible code. On the longer term, an extension of our method to include long 
wavelength acoustics would extend the validity of this method to M ~ 1 (see, for example, 
Gatti-Bono fc Colellall2006l ). We will also work on incorporating a flame model to capture 
the ignition of the flame (s) on the grid, so that we may continue the convective calculation 
in the presence of these first flames. 
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Ax 




Fig. 1. — The Cartesian grid and spherical base state (shown here in 2-d for simplicity, 
using 2Ar = Ax). Here we represent the spherical base state as concentric shells (black 
curved lines). Since the base state is not aligned with the Cartesian grid, we need to map 
between the two configurations. The '+' symbols represent the Cartesian zone centers. In 
our mapping from the radial profile to the Cartesian grid, the zones marked with the 'x' 
symbol are assigned the value from the gray-shaded radial bin. 
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Fig. 2. — The initial model used for the full star convection calculation. The top panel shows 
the density and the bottom panel shows the temperature. In both panels, the vertical dotted 
gray line represents the location of the low density cutoff — data outside of this cutoff are not 
used by our calculations. The vertical dashed gray line indicates where our sponge forcing 
term begins. For the temperature plot, the solid line represents the initial model used in 
our calculation and the dashed line represents the temperature structure for a completely 
isentropic model. 
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Fig. 3. — The inner sponge function, /damp, as a function of radius for p cuto fr = 3 x 10 6 g cm 
(solid line), and the outer sponge function for D = 5 x 10 8 cm and a 384 3 grid (dashed). 
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Fig. 4. — ($p)r/ (p)r vs - r f° r the test problem at 3 times. We see that the relative change in 
density resulting from a large amplitude velocity perturbation is small. 
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Fig. 5. — Development of the convective flow in the 384 3 calculation. Here we plot vorticity. 
The data scale is capped at 1.75 even though the maximum vorticity steadily climbs as 
the simulation progresses, reaching over 13 s~ l by the last panel. From left to right, top to 
bottom, the panels show the vorticity at 50, 100, 200, 400, 800, 1600, 3200, and 6400 s. 
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Fig. 6. — Radial velocity shown at 4 different times: a. 800 s; b. 3200 s; c. 3420 s; d. 7131.79 s. 
The latter time corresponds to the point of ignition. Red contours indicate outward moving 
fluid while blue contours indicate inward moving fluid. Two contour levels are used for each 
sign, ±1.2 x 10 6 cm s _1 and ±2 x 10 6 cm s _1 . The gray contour is a surface of constant 
density, p = p C utofr, marking the surface of the star. 



-33- 




Fig. 7. — The spherical angles, 9 and <p, computed from (v r )i as a function of time for the 
384 3 calculation. The vertical dotted lines represent where the angle <fi crosses the 0-360° 
boundary, and are not really discontinuities. We see that the direction of the dipole changes 
constantly throughout the simulation. 
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Fig. 8. — The peak radial velocity as a function of time for the 384 3 calculation. 
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Fig. 9. — The maximum temperature in the white dwarf as a function of time for the 384 3 
calculation. The temperature increase is highly nonlinear, ending at ignition. To show detail, 
we restrict the vertical range of the plot to 8 x 10 s K. The inset shows the structure of T pea k 
during the last ~200 s. We see large, but damped excursions in central temperature just 
prior to ignition. 




Fig. 10. — Average temperature as a function of radius in the white dwarf, shown at the 
initial time and 3 later times, for the 384 3 convection calculation. With time, the energy 
dumped into the star by reactions causes the temperature to increase throughout the star. 
The curve at 7132 s corresponds to the time of ignition. 
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Fig. 11. — ($p)r/(p)r vs. r for the 384 3 white dwarf convection problem at 3 times. The 
curve at 7132 s corresponds to the time of ignition. We see that at all times, {Sp) r /(p) r 
remains well below 1% everywhere inside the star. 



-38- 



8.0 



7.5 - 



7,0 - 



6.5 - 



T 



Pcaietl 
PaAuS 



:ixlfi 6 g cm -3 



[if 



g cm. 



i — r 





1000 2000 3000 4000 

t (s) 



5000 6000 7000 



Fig. 12. — The maximum temperature in the white dwarf as a function of time for the two 
different choices of p cuto fr (10 6 g cm -3 and 3 x 10 6 g cm" 3 ). Both simulations use a 256 3 grid. 
Here we see excellent agreement between the two cases, indicating that the peak temperature 
is insensitive to our choice of p cu toff- 
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Fig. 13. — The peak temperature, T pea k, as a function of time for the high resolution run 
(384 3 ) and two medium- resolution runs (256 3 ), offset so the time of ignition lines up. Here we 
show only the last 500 s leading up to ignition. We see that the temperature rise approaching 
ignition matches well between these different calculations. 
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Fig. 14. — The perturbational temperature (T — (T)) in three orthogonal slice planes (x-y, 
x-z, and y-z) passing through the point (2.48 x 10 8 cm, 2.49 x 10 8 cm, 2.52 x 10 8 cm), at 
a simulation time of 7131.79 s. Only the central 64 3 portion of the domain is shown (833.3 
km on a side). Ignition occurs in the zone where T — (T) is the largest. 
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Fig. 15. — The peak temperature, T pea k, in the white dwarf as a function of time for three 
different resolutions. We see that as we increase the resolution, the temperature increase is 
slower. The lowest resolution case reaches ignition very quickly. Once it ignites, the peak 
temperature climbs to ~ 10 10 K almost instantly. To show detail, we restrict the vertical 
range of the plot to 8 x 10 8 K. 
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Fig. 16. — The peak radial velocity, (i> r )peak> i n the white dwarf as a function of time for 
three different resolutions. 



